function [X, L, X_err, L_err, iter, flag] = ...
    cheby3(A, Xo, X_ref, L_ref, maxit, tol, print)

if (nargin==6)
    print = 1;
end

flag = 1;

% Reference values
% b   = -0.992856415680385/0.999999922699637;
% rho =  0.99403466159856/0.999999922699637;
%Xo = Xo/norm(Xo);
X0  = Xo;
X1	= Xo; 
X2  = Xo;
L2	= 1;

% cheb iterations using ref b and rho
p    = 0;
l    = 0;

Npi  = 5;  % # power iterations for updating rho
tNpi = 20;  % total PI count
Nci  = 20; % # chebyshev iterations per cycle 
tNci = 0;  % total CI count

[X2, L2, rho] = powermethodRHO(A, Xo, tNpi); 

R1 = ones(length(X1),1);
R2 = ones(length(X1),1);

iter = tNpi;
rho_index = 0;
while(iter <= maxit)
    
    if (p == 0)
        p = p + 1;
        iter = iter + Npi; 
        tNpi = tNpi + Npi;
        [X2, L2, rho] = powermethodRHO(A, X2, Npi);        
        b = -rho; % b is so close to rho in magnitude, this is probably fine
    end
    
    l = 0; % Chebyshev index
    while (l <= Nci && iter <= maxit)
        l    = l + 1;
        iter = iter + 1;  
        tNci = tNci + 1;
        
        % Save.
        X0 = X1;
        X1 = X2;
        R1 = R2;
        L1 = L2;
       
%         rho = 0.993972526759680;
%         b   = -0.992776504518967;
        % Get Chebyshev parameters.
        [alpha, beta] = alpha_beta(rho, b, l);
        
        % Update vectors and eigenvalue.
        V   = A * X1 ;   
        L2  = (V'*X1)/(X1'*X1); % Essentially a lagged Rayleigh
        X2  = X1 + alpha*(V/L2  - X0) + beta*(X1 - X0);
        %L2  = L1*norm(A*X2,2)/norm(A*X1,2);
        %X2 = X2 / norm(X2);
        % Internal convergence check.
        R2 = X2 - X1;
        %X2 = X2 / norm(X2);
       % q  = Q(R2, R1);
        q  = norm(R2)/norm(R1);
        r  = rho_est(rho, q, l);
        [restart,r_crit] = do_restart(q, rho, l);
        
        %disp(sprintf(' Restarting with rho = %6.4e (old rho = %6.4e)',r,rho));
        if (restart == 2)
            old_rho = rho;
           % rho = r;
%             if rho_index == 0
%                 rho = min([rho 0.9]);
%             elseif rho_index == 1
%                 rho = min([rho 0.95]);
%             elseif rho_index == 2
%                 rho = min([rho 0.985]);
%             end
%             if rho >= 1.0
%                % disp('warning: rho >= 1, setting to 0.99')
%                 rho = 0.99;
%             end
%            rho_index = rho_index + 1;
%            l   = 0;
            disp(sprintf(' Restarting with rho = %6.4e (old rho = %6.4e, l =%4i, r =%6.4e, log(q) =%6.4e)',r,old_rho,l,r_crit,q));
            break;
        end
        
        % Strictly for comparison studies
        %X     = X2/norm(X2,2); % ensure normed to 1
        
        % I've gotta be able to effect this.
        [X, L] = powermethod(A, X2, X_ref, L_ref, 2, tol, 0);
        
        X_err = norm(X - sign(X(1)/X_ref(1))*X_ref, 2);
        L_err = abs(L2 - L_ref);
       % printout(iter,X_err,L2,L1,print,rho)
        if ( X_err <= tol && L_err <= tol )
            flag = 0;
            break;
        end;
        if ( mod(iter, 20) == 0 )
            printout(iter,X_err,L2,L_err,print,rho)
        end
        
    end
    if (flag == 0)
        break;
    end
    p = 0; % let another power loop occur
end

if (flag==0 && print == 1)
    disp(' *** Chebyshev: Final Results *** ');
end
if (flag==1 && print==1)
    disp(' *** Warning: Chebyshev did not converge ***');
end
disp(sprintf(' Power Iters = %5i, Chebyshev Iters = %5i, TOTAL Iters = %5i', ...
    tNpi, tNci, iter));
printout(iter,X_err,L2,L_err,print,rho)
end

function printout(iter,X_err,L,L_err,print,rho)
if (print==1)
   disp(sprintf(...
       ' Iter = %5i, X_err = %6.5e, L = %6.5e, L_err = %6.5e, rho = %6.4e ',...
       iter, X_err,L,L_err,rho));
end
end

function [alpha, beta] = alpha_beta(rho, b, l)
% [alpha, beta] = alpha_beta(rho, b, l)
%
% Compute parameters for the Chebyshev recursion.  It is usually
% Here, b approximates lambda_N / lambda_1 and should be bounded from above
%  by that ratio.
% Here, rho approximates lambda_2 / lambda_1, the dominance ratio, and should be
%  bounded from below by that ratio.

% Reference values
%   b   = -0.992856415680385/0.999999922699637;
%   rho =  0.99403466159856/0.999999922699637;
if ( l == 1 )
    % special values for first iteration in cheby sequence
    alpha = 2 / (2 - rho - b);
    beta  = 0;
else
    y = (2.0 - rho - b) / (rho - b);
    alpha = 4.0 / (rho - b) * (Tn(l - 1, y) / Tn(l, y)); % Eqs. 11-2.12 in
    beta = Tn(l - 2, y) / Tn(l, y);                      % Hageman & Young
end

end

% Chebyshev polynomial of degree l evaluated at y
function z = Tn(l, y)
%if(abs(y) >= 1)
    z = cosh(l * acosh(y));
%else
%    z = cos(l * acos(y));
%end
end

% Estimate of rho
%   rho -- current estimate
%   Q   -- residual vector quotient
%   l   -- iteration index within current chebyshev cycle
function r = rho_est(rho, Q, l)
    if l == 1
        r = rho;
    else
        r = 0.5*rho * ( cosh( acosh(Q*Tn(l-1,2.0/rho-1.0)) / (l-1) ) + 1.0 );
    end
end

% Q-function from Hageman and Kellog.  Q is the residual vector
% quotient, and quantifies how much iterates are changing
function q = Q(R2,R1)
    q = norm(R2,2) / norm(R1,2);
end

% Restart check
function [y,R] = do_restart(Q, rho, l)
    y = 0; % Skip restart
    R = log(Q) / log( Tn(l-2,2/rho-1)/Tn(l-1,2/rho-1) );
    if ( R < 0.6 && l > 3 )
       y = 1; % Do restart if our convergence is sufficiently 
              % sub-optimal *and* we have generated a polynomial of at
              % least l
    end
end